function draw = drawnab(mu,sigma,a,b);
% function draw = drawnab(mu,sigma,a,b);
% Univariate sampling from a truncated normal distribution 
% in the interval [a,b]; 

if b <= a; 
    error('b must be greater than a'); 
end 
p1=normcdf(a,mu,sigma); 
p2=normcdf(b,mu,sigma); 
u=rand(1,1); 
te=norminv( u*(p2-p1) + p1 ); 
draw=mu + sigma*te ; 